using GeoPhyInv
using Statistics
using Plots; gr();
using ColorSchemesMedium
medium = Medium(:acou_homo2D, 5); # load a simple homogeneous acoustic medium from the gallery
update!(medium, [:vp, :rho], randn_perc = 5); # add some random noise to the medium
println(medium)> 2-D Medium with [:z, :x] in [[-1000.0, 1000.0], [-1000.0, 1000.0]]
> number of samples: [401, 401]
> sampling intervals: [5.0, 5.0]
> vp: min 1995.0397603607325 max 3076.4609700370515
> rho: min 1969.3995953763388 max 3034.211095525826
> Bounds:
5-element Named Vector{Vector{Float64}}
vcat │
──────┼───────────────────────────
:vp │ [2250.0, 2750.0]
:rho │ [2250.0, 2750.0]
:K │ [1.13906e10, 2.07969e10]
:KI │ [4.80841e-11, 8.77915e-11]
:rhoI │ [0.000363636, 0.000444444]AGeom
ageom = AGeom(medium.mgrid, :xwell, SSrcs(2)); # load a simple acquisition using `mgrid` of the medium
println(ageom)> 2-D acquisition with 2 supersource(s), [1, 1] source(s) and [10, 10] receiver(s)SrcWav
tgrid = range(0.0, stop = 2.0, length = 2500); # generate a time grid
wav = ricker(15.0, tgrid, tpeak = 0.25); # ricker wavelet
srcwav = SrcWav(tgrid, ageom, [:p]);
update!(srcwav, [:p], wav);Plotting the vp model
p1=heatmap(medium, :vp, seriescolor=cgrad(:roma))
scatter!(ageom, SSrcs())
scatter!(ageom, Recs())
plot(p1)Acoustic
Now we have all the required variables to create SeisForwExpt object and prepare the forward modeling. While creating, we switched the snaps_flag on, and instruct recording field at tsnaps. Once the Expt object is created, do the modeling without additional any memory allocations using update!
tsnaps=tgrid[1:div(length(tgrid),20):end] # store 20 snapshots
pa = SeisForwExpt(
FdtdAcoustic(),
medium = medium,
ageom = ageom,
srcwav = srcwav,
snaps_field = :p,
tsnaps = tsnaps,
tgrid = tgrid,
rfields = [:p],
verbose = true,
);[ Info: frequency bounds for propagating wavefield 1 are: [5.00e-01, 4.10e+01], with peak at: 1.05e+01
[ Info: spatial sampling (5.00e+00) can be as high as 5.49e+00
[ Info: time sampling (8.00e-04) can be as high as 1.29e-03Update pa to perform modelling; prints time and allocations
@time update!(pa);
FdtdAcoustic() supershot 1/2 1%|▏ | ETA: 0:02:18
FdtdAcoustic() supershot 1/2 23%|█████ | ETA: 0:00:20
FdtdAcoustic() supershot 1/2 29%|██████▍ | ETA: 0:00:18
FdtdAcoustic() supershot 1/2 34%|███████▍ | ETA: 0:00:16
FdtdAcoustic() supershot 1/2 38%|████████▌ | ETA: 0:00:15
FdtdAcoustic() supershot 1/2 43%|█████████▌ | ETA: 0:00:14
FdtdAcoustic() supershot 1/2 48%|██████████▋ | ETA: 0:00:12
FdtdAcoustic() supershot 1/2 53%|███████████▊ | ETA: 0:00:11
FdtdAcoustic() supershot 1/2 58%|████████████▉ | ETA: 0:00:09
FdtdAcoustic() supershot 1/2 64%|██████████████ | ETA: 0:00:08
FdtdAcoustic() supershot 1/2 69%|███████████████▏ | ETA: 0:00:07
FdtdAcoustic() supershot 1/2 74%|████████████████▍ | ETA: 0:00:06
FdtdAcoustic() supershot 1/2 79%|█████████████████▌ | ETA: 0:00:04
FdtdAcoustic() supershot 1/2 85%|██████████████████▋ | ETA: 0:00:03
FdtdAcoustic() supershot 1/2 90%|███████████████████▊ | ETA: 0:00:02
FdtdAcoustic() supershot 1/2 95%|█████████████████████ | ETA: 0:00:01
FdtdAcoustic() supershot 1/2 100%|██████████████████████| Time: 0:00:21
FdtdAcoustic() supershot 2/2 5%|█▏ | ETA: 0:00:18
FdtdAcoustic() supershot 2/2 11%|██▍ | ETA: 0:00:17
FdtdAcoustic() supershot 2/2 16%|███▍ | ETA: 0:00:17
FdtdAcoustic() supershot 2/2 20%|████▌ | ETA: 0:00:16
FdtdAcoustic() supershot 2/2 25%|█████▋ | ETA: 0:00:15
FdtdAcoustic() supershot 2/2 30%|██████▋ | ETA: 0:00:14
FdtdAcoustic() supershot 2/2 35%|███████▊ | ETA: 0:00:13
FdtdAcoustic() supershot 2/2 40%|████████▊ | ETA: 0:00:12
FdtdAcoustic() supershot 2/2 45%|█████████▊ | ETA: 0:00:11
FdtdAcoustic() supershot 2/2 49%|██████████▉ | ETA: 0:00:10
FdtdAcoustic() supershot 2/2 55%|████████████ | ETA: 0:00:09
FdtdAcoustic() supershot 2/2 60%|█████████████▏ | ETA: 0:00:08
FdtdAcoustic() supershot 2/2 65%|██████████████▍ | ETA: 0:00:07
FdtdAcoustic() supershot 2/2 71%|███████████████▌ | ETA: 0:00:06
FdtdAcoustic() supershot 2/2 76%|████████████████▊ | ETA: 0:00:05
FdtdAcoustic() supershot 2/2 81%|█████████████████▉ | ETA: 0:00:04
FdtdAcoustic() supershot 2/2 87%|███████████████████ | ETA: 0:00:03
FdtdAcoustic() supershot 2/2 92%|████████████████████▏ | ETA: 0:00:02
FdtdAcoustic() supershot 2/2 97%|█████████████████████▍| ETA: 0:00:01
FdtdAcoustic() supershot 2/2 100%|██████████████████████| Time: 0:00:19
timer output of methods inside fdtd time loop
──────────────────────────────────────────────────────────────────────────────
Time Allocations
────────────────────── ───────────────────────
Tot / % measured: 48.1s / 83.5% 720MiB / 12.0%
Section ncalls time %tot avg alloc %tot avg
──────────────────────────────────────────────────────────────────────────────
advance time step 5.00k 37.9s 94.4% 7.58ms 70.6MiB 81.9% 14.4KiB
record at receivers 5.00k 2.12s 5.27% 423μs 5.74MiB 6.67% 1.18KiB
add source 5.00k 132ms 0.33% 26.4μs 9.81MiB 11.4% 2.01KiB
──────────────────────────────────────────────────────────────────────────────
48.390481 seconds (12.70 M allocations: 761.451 MiB, 0.61% gc time, 18.07% compilation time)Extracting snaps
snaps1 = pa[:snaps, 1]; # extracting snaps of the first supersource
snaps2 = pa[:snaps, 2]; # second supersourcePlotting pressure snapshots after acoustic simulation
function plot_snapshots(name)
clip_perc = 90
pmax = maximum([maximum(pa[:snaps][it]) for it = 1:length(tsnaps)])
pmax -= pmax * 0.01 * clip_perc
anim = @animate for it = 1:length(tsnaps)
plot(
[
(
f = pa[:snaps, is][it];
heatmap(
f,
aspect_ratio = 1,
xlims = (1, size(f, 2)),
ylims = (1, size(f, 1)),
yflip = true,
seriescolor = cgrad(:seismic),
clims = (-pmax, pmax),
legend = false,
)
) for is = 1:2
]...,
title = string("snapshot at: ", round(tsnaps[it], digits = 5), " s"),
)
end
return gif(anim, string(name,".gif"), fps = 1)
end
plot_snapshots("acoustic_snaps")
Elastic
Similarly, we can do the elastic wave-equation modeling
medium = Medium(:elastic_homo2D, 5); # load a simple homogeneous elastic medium from the gallery
update!(medium, [:vp, :vs, :rho], randn_perc = 5); # add some random noise to the medium
println(medium)
ageom = AGeom(medium.mgrid, :xwell, SSrcs(2)); # load a simple acquisition using `mgrid` of the medium
println(ageom)> 2-D Medium with [:z, :x] in [[-1000.0, 1000.0], [-1000.0, 1000.0]]
> number of samples: [401, 401]
> sampling intervals: [5.0, 5.0]
> vp: min 1963.116052503446 max 3101.0698069092305
> vs: min 1154.7950463919065 max 1901.4136541003995
> rho: min 1949.1462081580023 max 3071.836020948651
> Bounds:
8-element Named Vector{Vector{Float64}}
vcat │
──────┼───────────────────────────
:vp │ [2250.0, 2750.0]
:vs │ [1350.0, 1650.0]
:rho │ [2250.0, 2750.0]
:K │ [5.92312e9, 1.08144e10]
:KI │ [9.24695e-11, 1.6883e-10]
:rhoI │ [0.000363636, 0.000444444]
:mu │ [4.10062e9, 7.48688e9]
:muI │ [1.33567e-10, 2.43865e-10]
> 2-D acquisition with 2 supersource(s), [1, 1] source(s) and [10, 10] receiver(s)Add source term on :vz grid
srcwav = SrcWav(tgrid, ageom, [:vz]);
update!(srcwav, [:vz], wav);Perform modelling
pa = SeisForwExpt(
FdtdElastic(),
medium = medium,
ageom = ageom,
srcwav = srcwav,
snaps_field = :vz,
tsnaps = tsnaps,
rfields = [:vz],
tgrid = tgrid,
verbose = true,
)
@time update!(pa);
snaps = pa[:snaps, 1];
snaps = pa[:snaps, 2];[ Info: frequency bounds for propagating wavefield 1 are: [5.00e-01, 4.10e+01], with peak at: 1.05e+01
┌ Warning: decrease maximum spatial sampling (5.00e+00) below 2.82e+00
└ @ GeoPhyInv ~/work/GeoPhyInv.jl/GeoPhyInv.jl/src/fdtd/stability.jl:38
[ Info: time sampling (8.00e-04) can be as high as 1.10e-03
FdtdElastic() supershot 1/2 2%|▍ | ETA: 0:01:02
FdtdElastic() supershot 1/2 4%|█ | ETA: 0:00:46
FdtdElastic() supershot 1/2 7%|█▌ | ETA: 0:00:42
FdtdElastic() supershot 1/2 9%|██▏ | ETA: 0:00:40
FdtdElastic() supershot 1/2 12%|██▋ | ETA: 0:00:39
FdtdElastic() supershot 1/2 14%|███▎ | ETA: 0:00:37
FdtdElastic() supershot 1/2 16%|███▊ | ETA: 0:00:36
FdtdElastic() supershot 1/2 19%|████▍ | ETA: 0:00:35
FdtdElastic() supershot 1/2 21%|████▉ | ETA: 0:00:34
FdtdElastic() supershot 1/2 24%|█████▍ | ETA: 0:00:33
FdtdElastic() supershot 1/2 26%|█████▉ | ETA: 0:00:32
FdtdElastic() supershot 1/2 28%|██████▌ | ETA: 0:00:31
FdtdElastic() supershot 1/2 30%|███████ | ETA: 0:00:30
FdtdElastic() supershot 1/2 33%|███████▌ | ETA: 0:00:30
FdtdElastic() supershot 1/2 35%|████████ | ETA: 0:00:29
FdtdElastic() supershot 1/2 37%|████████▌ | ETA: 0:00:28
FdtdElastic() supershot 1/2 39%|█████████ | ETA: 0:00:27
FdtdElastic() supershot 1/2 41%|█████████▌ | ETA: 0:00:26
FdtdElastic() supershot 1/2 43%|██████████ | ETA: 0:00:25
FdtdElastic() supershot 1/2 46%|██████████▌ | ETA: 0:00:24
FdtdElastic() supershot 1/2 48%|███████████▏ | ETA: 0:00:23
FdtdElastic() supershot 1/2 51%|███████████▋ | ETA: 0:00:22
FdtdElastic() supershot 1/2 53%|████████████▏ | ETA: 0:00:21
FdtdElastic() supershot 1/2 55%|████████████▊ | ETA: 0:00:20
FdtdElastic() supershot 1/2 58%|█████████████▍ | ETA: 0:00:18
FdtdElastic() supershot 1/2 61%|█████████████▉ | ETA: 0:00:17
FdtdElastic() supershot 1/2 63%|██████████████▌ | ETA: 0:00:16
FdtdElastic() supershot 1/2 66%|███████████████▏ | ETA: 0:00:15
FdtdElastic() supershot 1/2 68%|███████████████▋ | ETA: 0:00:14
FdtdElastic() supershot 1/2 71%|████████████████▎ | ETA: 0:00:13
FdtdElastic() supershot 1/2 73%|████████████████▉ | ETA: 0:00:12
FdtdElastic() supershot 1/2 76%|█████████████████▍ | ETA: 0:00:10
FdtdElastic() supershot 1/2 78%|██████████████████ | ETA: 0:00:09
FdtdElastic() supershot 1/2 81%|██████████████████▋ | ETA: 0:00:08
FdtdElastic() supershot 1/2 83%|███████████████████▏ | ETA: 0:00:07
FdtdElastic() supershot 1/2 86%|███████████████████▊ | ETA: 0:00:06
FdtdElastic() supershot 1/2 88%|████████████████████▍ | ETA: 0:00:05
FdtdElastic() supershot 1/2 91%|████████████████████▉ | ETA: 0:00:04
FdtdElastic() supershot 1/2 94%|█████████████████████▌ | ETA: 0:00:03
FdtdElastic() supershot 1/2 96%|██████████████████████▏| ETA: 0:00:02
FdtdElastic() supershot 1/2 99%|██████████████████████▊| ETA: 0:00:01
FdtdElastic() supershot 1/2 100%|███████████████████████| Time: 0:00:42
FdtdElastic() supershot 2/2 3%|▋ | ETA: 0:00:38
FdtdElastic() supershot 2/2 5%|█▏ | ETA: 0:00:38
FdtdElastic() supershot 2/2 8%|█▊ | ETA: 0:00:37
FdtdElastic() supershot 2/2 10%|██▎ | ETA: 0:00:37
FdtdElastic() supershot 2/2 12%|██▉ | ETA: 0:00:36
FdtdElastic() supershot 2/2 15%|███▍ | ETA: 0:00:35
FdtdElastic() supershot 2/2 17%|████ | ETA: 0:00:34
FdtdElastic() supershot 2/2 19%|████▌ | ETA: 0:00:34
FdtdElastic() supershot 2/2 22%|█████ | ETA: 0:00:33
FdtdElastic() supershot 2/2 24%|█████▌ | ETA: 0:00:32
FdtdElastic() supershot 2/2 26%|██████▏ | ETA: 0:00:31
FdtdElastic() supershot 2/2 29%|██████▋ | ETA: 0:00:30
FdtdElastic() supershot 2/2 31%|███████▏ | ETA: 0:00:29
FdtdElastic() supershot 2/2 33%|███████▋ | ETA: 0:00:28
FdtdElastic() supershot 2/2 36%|████████▏ | ETA: 0:00:28
FdtdElastic() supershot 2/2 38%|████████▋ | ETA: 0:00:27
FdtdElastic() supershot 2/2 40%|█████████▏ | ETA: 0:00:26
FdtdElastic() supershot 2/2 42%|█████████▊ | ETA: 0:00:25
FdtdElastic() supershot 2/2 45%|██████████▎ | ETA: 0:00:24
FdtdElastic() supershot 2/2 47%|██████████▊ | ETA: 0:00:23
FdtdElastic() supershot 2/2 49%|███████████▍ | ETA: 0:00:22
FdtdElastic() supershot 2/2 52%|███████████▉ | ETA: 0:00:21
FdtdElastic() supershot 2/2 54%|████████████▌ | ETA: 0:00:20
FdtdElastic() supershot 2/2 57%|█████████████▏ | ETA: 0:00:18
FdtdElastic() supershot 2/2 60%|█████████████▊ | ETA: 0:00:17
FdtdElastic() supershot 2/2 62%|██████████████▎ | ETA: 0:00:16
FdtdElastic() supershot 2/2 65%|██████████████▉ | ETA: 0:00:15
FdtdElastic() supershot 2/2 67%|███████████████▌ | ETA: 0:00:14
FdtdElastic() supershot 2/2 70%|████████████████▏ | ETA: 0:00:13
FdtdElastic() supershot 2/2 73%|████████████████▊ | ETA: 0:00:12
FdtdElastic() supershot 2/2 75%|█████████████████▎ | ETA: 0:00:10
FdtdElastic() supershot 2/2 78%|█████████████████▉ | ETA: 0:00:09
FdtdElastic() supershot 2/2 80%|██████████████████▌ | ETA: 0:00:08
FdtdElastic() supershot 2/2 83%|███████████████████▏ | ETA: 0:00:07
FdtdElastic() supershot 2/2 86%|███████████████████▋ | ETA: 0:00:06
FdtdElastic() supershot 2/2 88%|████████████████████▎ | ETA: 0:00:05
FdtdElastic() supershot 2/2 91%|████████████████████▉ | ETA: 0:00:04
FdtdElastic() supershot 2/2 93%|█████████████████████▌ | ETA: 0:00:03
FdtdElastic() supershot 2/2 96%|██████████████████████ | ETA: 0:00:02
FdtdElastic() supershot 2/2 99%|██████████████████████▋| ETA: 0:00:01
FdtdElastic() supershot 2/2 100%|███████████████████████| Time: 0:00:41
timer output of methods inside fdtd time loop
──────────────────────────────────────────────────────────────────────────────
Time Allocations
────────────────────── ───────────────────────
Tot / % measured: 85.8s / 96.8% 233MiB / 47.6%
Section ncalls time %tot avg alloc %tot avg
──────────────────────────────────────────────────────────────────────────────
advance time step 5.00k 80.9s 97.4% 16.2ms 109MiB 98.4% 22.3KiB
record at receivers 5.00k 2.09s 2.52% 419μs 922KiB 0.81% 189B
add source 5.00k 50.7ms 0.06% 10.1μs 922KiB 0.81% 189B
──────────────────────────────────────────────────────────────────────────────
85.972455 seconds (4.13 M allocations: 251.794 MiB, 0.16% gc time, 3.52% compilation time)Plotting is vs model
p1=heatmap(medium, :vs, seriescolor=cgrad(:roma))
scatter!(ageom, SSrcs())
scatter!(ageom, Recs())
plot(p1)Plotting snapshots after elastic simulation
plot_snapshots("elastic_snaps") # P and S waves!